Microscopic origin of the structural phase transitions at the Cr 2 Os (0001) surface 
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The surface of a Cr2C>3 (0001) film epitaxially grown on Cr undergoes an unusual reentrant se- 
quence of structural phase transitions (lxl-4 y/3 x \/3 -+lxl). In order to understand the 
underlying microscopic mechanisms, the structural and magnetic properties of the Cr2 03 (0001) sur- 
face are here studied using first-principles electronic structure calculations. Two competing surface 
Cr sites are identified. The energetics of the surface is described by a configurational Hamiltonian 
with parameters determined using total energy calculations for several surface supercells. Effects of 
epitaxial strain and magnetic ordering on configurational interaction are also included. The ther- 
modynamics of the system is studied using Monte Carlo simulations. At zero strain the surface 
undergoes a 1 x 1 — + v3 x \/3 ordering phase transition at Tq ~ 165-RT. Tensile epitaxial strain 
together with antiferromagnetic ordering drive the system toward strong configurational frustration, 
suggesting the mechanism for the disordering phase transition at lower temperatures. 
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I. INTRODUCTION 

Metal oxides demonstrate a variety of physical and 
chemical properties, sometimes in intriguing combina- 
tions. Apart from being ubiquitous in nature, metal- 
oxide surfaces and interfaces find diverse technological 
applications and are being explored for potential use 
in future electronic devices. In particular, surfaces of 
magnetoelectric antiferromagnets such as &2O3 possess 
an equilibrium surface magnetization^— making them 
suitable for use as active layers in electrically-switchable 
magnetic nanostructuresji 

The Cr 2 03 (0001) surface has been a subject of many 
experimental 5 - - — and theoretical 5 ^^ - — studies, but its 
surface remains poorly understood. Low-energy electron 
diffraction (LEED) experiments for a thin &2O3 (0001) 
film grown on a Cr (110) single crystal revealed an un- 
usual reentrant structural phase transition, 5 in which the 
surface structure changes from 1 x 1 to vo x yo and back 
to 1 x 1 under cooling from room temperature to 150 K 
and then further down to 100 K. The origin of these 
phase transitions is not understood. While the high- 
temperature transition may, as suggested by the LEED 
data* 5 - be a conventional order-disorder transition, the 
second one is unusual in that a more symmetric phase 
appears at lower temperatures. 

The situation is further complicated by the fact that, 
as shown by Takano et air , both phase transitions disap- 
pear for thicker Cr 2 03 films grown in a similar way. This 
suggests that the epitaxial strain has an important effect 
on the surface energetics. Oxidation of the Cr 110 surface 
was investigated by LEED and Auger spectroscopy,— 
and it was found that under growth conditions similar 
to those of Ref. 5 the thin Cr 2 3 (0001) film is subject 
to a tensile epitaxial strain of about 1.5%. 

In this paper we study the structure of the &2O3 
(0001) surface using first-principles electronic structure 
calculations and Monte Carlo simulations. Our results 
suggest that the dynamics of the system is driven by 
the occupation of two competing surface Cr sites. The 



system can be mapped to an Ising model on a two- 
dimensional hexagonal lattice in external field. The first 
phase transition is clearly identified as a conventional or- 
dering transition, and the theoretical transition temper- 
ature is found to be in good agreement with experiment. 
Our calculations further reveal a strong effect of tensile 
epitaxial strain, which parametrically drives the system 
towards configurational frustration, particularly in com- 
bination with antiferromagnetic ordering. An explana- 
tion of the second phase transition is offered based on 
these results. 

The paper is organized as follows. In Section [TT] we de- 
scribe the computational methods. Section Mil presents 
the results on the configurational and magnetic energetics 
of the Cr203 (0001) surface, including the identification 
of the competing surface Cr sites, the construction of the 
configurational Hamiltonian, the analysis of magnetic in- 
teractions, and the evaluation of the ground-state phase 
diagram. Section IIVI deals with configurational thermo- 
dynamics of the surface, and Section [V] discusses the re- 
lation of the results to experiments. The electronic struc- 
ture of the Cr203 (0001) surface is presented in Section 
IVI1 and its magnetic properties in Section IVII1 Section 
IVIIII draws the conclusions. 



II. COMPUTATIONAL METHODS 

Electronic structure calculations were performed using 
the projector-augmented wave method^ implemented in 
the VASP code-i^ For the Cr 3d shell we employed the 
rotationally-invariant LSDA+J7 method^ with (7 = 4 
eV and J = 0.58 eV. This method was preferred over 
GGA+U adopted in Ref. [l6| due to its better descrip- 
tion of the structural, electronic, and magnetic properties 
of bulk Cr 2 03^ Different surface superstructures were 
modeled using supercells representing symmetric slabs 
with eight atomic layers of O and 16 atomic layers of Cr 
stacked along the (0001) direction. The periodically re- 
peating slab is separated from its image by 1.5 nm of vac- 
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uum. We considered lxl, 1x2, 1x3, and y/3x y3 surface 
supercells (where lxl corresponds to the hexagonal unit 
cell of bulk C^Oa). The lateral dimensions of the un- 
strained supercell were fixed to the calculated equilibrium 
bulk values^ for the strained case these values were used 
as a reference. Apart from these constraints, the ionic po- 
sitions were relaxed until the Hellmann-Feynman forces 
were converged to less than 0.01 eV/A. The plane-wave 
energy cutoff was fixed to 520 eV and the Brillouin zone 
integration was performed using T-centered Monkhorst- 
Pack grids^ For relaxation we used Gaussian smearing 
of 0.1 eV and a /c-point mesh equivalent to or denser than 
4 x 4 x 1 for the lxl surface supercell. We checked the 
convergence with respect to the number of Appoints, the 
energy cutoff for the plane wave expansion, the size of the 
vacuum region, and the thickness of the slab. These tests 
indicate that the total energies are generally converged 
to within 1 meV. Density of states (DOS) calculations 
were performed using Gaussian smearing of 0.02 eV and 
a fc-point mesh equivalent to or denser than 8 x 8 x 1 for 
the lxl supercell. 

The energy barriers for the thermally-activated jump- 
ing of Cr ions between the two competing surface sites 
were calculated using the nudged elastic band method.— 
Seven images were inserted between the two energy min- 
ima, and in each image the ions were relaxed so that 
forces perpendicular to the reaction path were smaller 
than 0.05 eV/A. 



III. SURFACE ENERGETICS 

Cr203 crystalizes in the corundum structure with the 
R3c space group. It can be viewed as a stacking of buck- 
led honeycomb Cr double layers along the (0001) direc- 
tion with quasi- hexagonal closed-packed O layers in be- 
tween, see Fig.[TJ The (0001) surface is polar, and simple 
electrostatic arguments suggest that non-stoichiometric 
terminations by an O layer or by a Cr double layer should 
lead to divergent electrostatic potential in the bulk. On 
the other hand, the surface can terminate in the mid- 
dle of the buckled Cr layer so that only half of the Cr 
ions from this layer remain on the surface. Although still 
polar, this termination is stoichiometric, and the electro- 
static potential in the bulk is not divergent. It can there- 
fore be expected that this termination is energetically 
favorable. Indeed, surface termination by a single Cr 
layer was consistent with LEED 6 and scanning tunneling 
microscope^ measurements of the Cr203 (0001) surface 
in ultrahigh vacuum. Further, first principles calcula- 
tions by Rohrbach et al*& based on the GGA+C/ method 
have shown that this termination has the lowest surface 
energy (compared to all others considered) over the entire 
range of oxygen chemical potential where Cr203 is sta- 
ble. Note that earlier results based on the GGA method, 
which leads to grossly incorrect electronic and magnetic 
properties^ were quite different, I 5 In this work we only 
consider the Cr203 (0001) surface terminated by a single 



layer of Cr. 

A. Surface sites 

The location of the Cr ions within the single Cr ter- 
minating layer has been debated. Within the double 
Cr layer there are three possible octahedral sites, two 
of them being occupied in the bulk. They give rise to 
three nonequivalent surface sites (A, C and D, see Fig. 
[IJ that surface Cr ions can occupy. Occupation of site 
A corresponds to the continuation of the bulk structure. 
Further, as pointed out by Gloege et. air- the surface Cr 
ion can jump below the oxygen subsurface layer and oc- 
cupy the empty octahedral site within the underlying Cr 
double layer. This interstitial site is directly underneath 
the surface site A, and we denote it by B (see Fig. [1]). 




FIG. 1. Slab geometries for four considered lxl surface ter- 
minations. Gray and red spheres represent Cr and O atoms, 
respectively. 

In order to identify the energetically preferable sites, 
we therefore considered four lxl surface models cor- 
responding to the exclusive occupation of sites A, B, C, 
or D, respectively. In all cases a significant inward re- 
laxation was observed, as expected for a nominally polar 
surface. The relaxation data for models A and B are 
included in Appendix [A] We define the surface energy as 

E s = — ( Egiab — r^—Ebuik) /N s (1) 

2 V i^bulk J 

Here E s i a b is the ground state energy of the slab for the 
given surface model with magnetic structure correspond- 
ing to bulk Cr203, E^uik is the ground state energy per 
unit cell of bulk Cr203, N s i a b and N^uik are the num- 
bers of atoms in the slab and in the bulk unit cell, and 
N s is the number of surface Cr atoms on one side of the 
slab. The surface energies for the four lxl surface ter- 
minations are given in Table [I] The surface energy is the 
lowest when site A is occupied. Occupation of sites C and 
D leads to much higher surface energies, and we therefore 
do not consider their occupation in the subsequent anal- 
ysis. On the other hand, the surface energy of model B is 
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TABLE I. Surface energies of 1 x 1 surface models with dif- 
ferent surface sites occupied. 





A 


B 


C 


D 


E„, eV 


2.909 


3.077 


5.002 


5.847 



only slightly higher than that of model A. Thus, sites A 
and B can both be partially occupied, which can lead to 
non-trivial ordered terminations and phase transitions; 
these issues are addressed in the following subsections. 

The identification of sites A and B as the most fa- 
vorable agrees with LEED measurements and molecular 
dynamics simulations of Ref. [gj, as well as with surface 
X-ray diffraction data^ but recent LEEDiii and surface 
X-ray diffraction (SXRD)ii studies have questioned the 
single Cr layer surface termination and reached different 
conclusions. In Ref. 1Q a non-stoichiometric surface with 
a partial occupation of four Cr layers near the surface 
was obtained, but the best-fit i?-factor R p = 0.48 was 
poor, as noted by the authors. In Ref. [ill the best fit for 
the SXRD measurements was obtained for a surface ter- 
minated with a partially occupied double Cr layer (sites 
A and C) and one more partially occupied Cr layer below 
that. Partial occupancy of site C is difficult to reconcile 
with the very high surface energy of surface model C (2 
eV per Cr site higher compared to model A), although 
this site could, in principle, be stabilized by intersite in- 
teractions or by depletion of Cr atoms in the subsurface 
Cr layers. Since the configurational models required to 
explore such unconventional terminations would be very 
complicated, we did not attempt to consider them. We 
also note that the site occupations must be integer in the 
ground state. Further analysis may be required as more 
experimental evidence becomes available. 

While the spin-orbit coupling in bulk Cr 2 03 is small^ 
the reduced coordination could make it more important 
at the surface. To estimate its role, we calculated the 
energy difference between surface models A and B in the 
presence of spin-orbit coupling (taking the structures re- 
laxed without it). It was found that spin-orbit coupling 
changes this energy difference by 0.4 meV. This energy is 
small compared to all important structural and exchange 
interaction parameters, and therefore spin-orbit coupling 
was neglected in all subsequent calculations. 



B. Configurational interaction at the surface 

Surface A sites form a two-dimensional hexagonal lat- 
tice, and there is a B site directly underneath every A 
site. Based on the surface energies calculated in the pre- 
vious subsection, we assume that at every hexagonal lat- 
tice site the Cr atom occupies either site A or site B. 
Therefore, we can introduce an occupation number rij, 
where i denotes a 2D hexagonal lattice site, such that rii 
is equal to 1 if site B is occupied and if site A is oc- 
cupied. The following configurational Hamiltonian can 



therefore be introduced: 

n = V int ({m}) + hJ2ni (2) 

i 

The first term includes the configurational interaction 
between surface Cr ions, and the second term takes into 
account that sites A and B are inequivalent. Since the 
total number of A sites is not conserved, this Hamilto- 
nian is isomorphic to an interacting Ising model on a 2D 
hexagonal lattice in external magnetic field. 

The introduction of the 2D hexagonal lattice is based 
on the spatial arrangement of A sites on the surface. 
Note, however, that the true symmetry of the (disor- 
dered) Cr2C>3 surface is lower: apart from the transla- 
tions, there are only C3 axes passing through the Cr sites. 
To take this difference into account, one can formally as- 
sign a direction to each bond on the 2D hexagonal lattice. 
The directions of the six nearest-neighbor bonds should 
be made alternating (i. e. three incoming and three out- 
going bonds). The directionality of the bonds can be 
reflected in the interaction term in the Hamiltonian 
For example, the pair interaction parameter may be dif- 
ferent for a bond pointing from site A to site B and for 
a bond pointing from site B to site A. However, we are 
mainly interested in the total energies of different con- 
figurations {rii} which are not strongly affected by the 
directionality of the bonds. In fact, it can be shown that 
for pairwise interaction of any range the total energies 
do not depend on whether the bond directionality is in- 
cluded or not. (This is because the total numbers of 
A— >B and B— >A bonds are always equal in all coordi- 
nation spheres.) Even when many-body interactions are 
present, the contribution of the bond directionally to the 
total energy is zero for most ordered configurations. In 
particular, among all configurations shown in Fig. [5] only 
A7B5 (6 x 6) has a nonzero contribution, but thus struc- 
ture is not important for any of the following. Moreover, 
as we will show below, the surface structure is governed 
by non-directional electrostatic interactions. We there- 
fore do not introduce any non-directional terms in the 
Hamiltonian, which makes the assignment of bond direc- 
tions superfluous. 

In order to proceed, we use the cluster expansion 
approach^ - — which is widely used in the studies of 
bulk alloy thermodynamics. Specifically, we need to 
adopt some particular representation of Vi n t({rii}) and 
fit it to the calculated total energies of different ordered 
configurations {n^}. However, due to large size of the 
system the calculations are only feasible for a few rela- 
tively small supercells (see Table|TT|), and we must request 
that Vi n t{{rii\) has but a small number of parameters. 
We construct such a representation based on physical 
grounds (rather than trial-and-error) and then validate 
the results by the quality of the fit. 

The structure of polar surfaces is expected to be dom- 
inated by electrostatic interactions. For example, for the 
polar GaAs (001) surface it was shown that surface en- 
ergy differences between different orderings are well de- 
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scribed by a simple electrostatic model^ 9 . We therefore 
include electrostatic interaction in Vi nt ({rii}) by treating 
surface Cr ions as point charges q interacting via classi- 
cal Coulomb forces screened by a dielectric constant e. 
We started by assuming that the positions of sites A and 
B do not depend on the environment, but this simple 
model was found to be inaccurate. However, it can be 
significantly improved by including the effect of atomic 
relaxations. 

For a given configuration {rij} the total energy can be 
reduced by shifting of the surface Cr ions from their aver- 
age positions at sites A and B. Such relaxation terms are 
often important in the thermodynamics of strongly size- 
mismatched bulk alloys. Although the Cr lattice sites 
at the Cr203 surface are located rather far from each 
other, we found that the vertical (normal to the surface) 
coordinate of a surface Cr ion occupying site A depends 
rather strongly on its environment (the shift can be as 
large as 0.4 A, see Appendix |A"|) . Other ions, including 
Cr atoms at site B, shift much less, and we therefore only 
consider relaxations of the A sites. (This approximation 
is justified by the resulting high quality of the fitting.) 
We introduce a vertical coordinate z, for each occupied 
A site and minimize it for the given configuration {rii}. 
(Thus, Zi are treated as adiabatically "fast" variables.) 
The electrostatic interaction contributes a vertical force 
depending on the occupation numbers at other sites of 
the lattice. The contribution to the total energy depend- 
ing on Zi is written as 

H(n, z) = -7 n% Oi - z f - ~ ^ 'zjf.^i ( 3 ) 

i ij l 3 

Here we defined hi = 1 — n*. The first term represents the 
elastic contribution for each A site as a simple harmonic 
oscillator with stiffness 7 and equilibrium position zq . 
The second term describes the electrostatic interaction 
with B sites. (Small vertical forces from other A sites 
are neglected.) Here dij is the distance between sites i 
and j, and the B sites are assumed to lie at z = 0. Since 
Zi <C dij, we have used the dipole approximation with 

Pi = qzi. 

Introducing a small parameter a = q 2 /(jea 3 ) 1 where 
a is the 2D hexagonal lattice parameter, and minimizing 
(J3]) with respect to z,-, we obtain, to first order in a: 

Zi = z + az Ri (4) 

where Ri = X^j^jCij an d Cij = (a/dij) 3 - This relation 
agrees perfectly with relaxation data for surface super- 
cells (see Appendix [A"| . Substituting z, in ([3]), we obtain 
to first order in a: 

V* ni ({n,}) = -^V^Cijnihj - X^^R 2 (5) 

ij i 

where V = q 2 z 2 /(ea 3 ) and X = aV/2. The first two- 
body term represents the dipolar interactions assuming 
fixed A site positions Zi = zq. The second three-body 



term is the lowest-order correction due to the A site 
shifts. Note that the parameter V is positive, because 
an unlike AB bond is longer than an AA or a BB bond 
due to the vertical shift, and all Cr ions are positively 
charged. This transparent physical mechanism generates 
an ordering tendency in our system. 

The resulting configurational Hamiltonian contains 
three parameters: h, V and X. Note that Vint vanishes 
when all rij = or all n,; = 1, and therefore h gives the 
positive energy difference between models A and B in 
Table H 

As noted in the Introduction, thin films of &2O3 
demonstrating phase transitions are subject to a tensile 
epitaxial strain of about 1.5%. Therefore, in the follow- 
ing we consider two cases: (1) unstrained surface, and 
(2) surface subject to a 1.5% in-plane tensile strain. 

The three parameters of the configurational model are 
fitted to the calculated surface energies of several ordered 
configurations listed in Table [TT1 and illustrated in Fig. [2] 
The standard take-one-out cross-validation (CV) scored 
is used to evaluate the predictive power of the fit. It can 
be seen that the model ([5]) provides an excellent fitting to 
the calculated energies for both unstrained and strained 
surfaces, supporting our assumptions about the physical 
interaction mechanisms. Note that the parameter X is 
small compared to V in agreement with our assumptions. 
Nevertheless, the three-body term is essential for obtain- 
ing a good fit (see Appendix IB1 for further discussion). 
Without this term a much larger CV score is obtained, 
and the take-one-out prediction for model B is about 50 
meV off. The quality of the fit is also significantly im- 
paired if the range of the electrostatic interaction is cut 
off in real space. 

One might expect that Vi nt ({rii}) for nearest neighbors 
could also have a contribution of non-electrostatic origin. 
However, since the surface Cr ions are rather far from 
each other (a w 5 A) and the surface remains insulating, 
this contribution should be short-ranged and relatively 
small. Indeed, the addition of a nearest-neighbor pair 
or three-body (triangle) interaction to Vi n t({n.i}) did not 
improve the quality of the fit. 

The main effect of tensile strain on the configurational 
Hamiltonian is the decrease of the parameter h by about 
a factor of two compared to the unstrained surface. This 
effect can be understood by noting that h represents the 
local preference of the bulk-like surface site A over the in- 
terstitial site B. Under tensile strain the lattice expands, 
leaving more space available for the Cr stom at site B. 
This reduces the interstitial pressure and thereby the en- 
ergy cost of occupying site B. 



C. Effect of magnetic ordering 

So far we have discussed the fitting of the configura- 
tional Hamiltonian to surface energies for antiferromag- 
netically ordered supercells. The directions of the local 
moments at B sites were assigned similar to A sites, con- 
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TABLE II. Ground state (AFM) and paramagnetic (PM) sur- 
face energies for different configurations for the Cr203 (OOOf ) 
surface under zero strain and under 1.5% tensile epitaxial 
strain. A x ~B y (di x da) denotes the configuration with a sur- 
face supercell spanned by vectors of length d\ and d 2 and 
containing x (y) surface Cr ions in position A (B). The con- 
figurations in the table are shown explicitly in Fig. [3] The 
values are given with respect to the surface energy of model 
A. Corresponding fitted values of parameters of the configu- 
rational Hamiltonian together with the misfit and the average 
cross-validation score are also given. Further, the critical tem- 
perature of the (\/3 x \/3) to (1 x 1) order-disorder transition 
obtained from MC is also given. The surface energies, param- 
eters of the configurational Hamiltonian, the misfit and the 
average cross-validation score are all given in meV while the 
critical temperature is in K. 





Unstrained 


Strained 


AFM 


PM 


AFM 


PM 


B (1 x 1) 


168 


160 


77 


91 


AB (1 x 2) 


-51 


-48 


-71 


-64 


A 2 B (1 x 3) 


-42 


-41 


-55 


-51 


AB 2 (1 x 3) 


-2 


-2 


-41 


-31 


A 2 B (V3 x v/3) 


-62 


-57 


-70 


-65 


AB 2 (VS x v/3) 


-28 


-22 


-60 


-49 


h 


168 


160 


76 


90 


V 


64 


64 


52 


56 


X 


1.6 


1.4 


1.2 


1.0 


misfit 


1 


1 


1 


1 


CV 


2 


7 


3 


3 


T C 


165 ±5 


165 ±5 




50 ± 10 



tinuing the bulk antiferromagnetic structure. The surface 
magnetic structure may, however, be different from the 
bulk one. We checked this by recalculating the surface 
energies for different magnetic configurations of a few Cr 
sites closest to the surface. These sites included A and B 
sites, as well as the two Cr sites in the underlying buckled 
honeycomb Cr layer (types 2 and 3 in the order of depth, 
see Fig. [3]). Assuming that Cr ions of the same type 
always have the same spin direction, for each input sur- 
face configuration we calculated the total energy for all 
possible configurations of the four near-surface Cr sites 
(A, B, 2, and 3), while keeping the rest of the slab in its 
bulk magnetic structure. 

We found that the lowest surface energy corresponds 
to the continuation of the bulk magnetic structure for 
all surface configurations with the exception of surface 
model B (1 x 1). For this model the surface energy is re- 
duced by 43 meV by flipping of the local moment on site 
2, thereby making it parallel to those on sites B and 3. 
Nevertheless, the surface energy of model B listed in the 
first column of Table U and used in the fitting corresponds 
to the continuation of the bulk structure. This preserves 
consistency with the surface energies of other configu- 
rations. The effect of this choice on thermodynamics is 
small, because surface model B has a large surface energy. 
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FIG. 2. Surface configurations mentioned in the text. Config- 
urations (a)-(g) were used to fit the parameters of the configu- 
rational Hamiltonian. Configurations (a), (c), (f) and (h)-(n) 
have been identified as possible ground states. 



Magnetic disorder present at finite temperatures may 
affect the relative energies of different surface models and 
thereby influence the thermodynamic properties. A com- 
plete solution requires that the structural (m) and mag- 
netic degrees of freedom are both included in the effec- 
tive Hamiltonian. We did not attempt to construct such 
a Hamiltonian, but rather considered the effect of com- 
plete magnetic disorder in the paramagnetic phase on the 
structural interaction parameters. To this end, for each 
of the input surface models we have fitted the surface 
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energies to a surface Heisenberg Hamiltonian 

n = - \ J2 ■ *i - E H ^ + N ° E * M ( 6 ) 

ij i 

Here summation runs over Cr ions belonging to one of 
the four types defined above (A, B, 2, 3), and Sj is a 
unit vector parallel to the local moment of the i-th ion. 
Assuming that the exchange coupling does not extend 
further than in the bulk, the only nonzero exchange pa- 
rameters are those between nearest-neighbor Cr ions of 
different types (see Fig. [3] for an illustration) . In addi- 
tion, each Cr ion interacts with an effective exchange field 
Hi set up by the bulk. We assumed that Cr ions of the 
same type are equivalent. Under these assumptions the 
number of parameters reduces to 10: Ha, Hb, H 2 , H 3 , 
Ja2, J A3, Jb2, Jb3, J23, and Ef M . For lxl surface 
models only seven of these parameters remain. We have 
fitted these parameters using 8 magnetic configurations 
for the lxl surface models and using 16 magnetic con- 
figurations for models with a larger unit cell. Good fits 
were obtained for all surface models. Table Hill lists the 
fitted parameters and the misfits for surface models A, 
B, and A2B. These parameters will be discussed further 
in Section Ivnl 




FIG. 3. Magnetic model for the (0001) Cr 2 3 surface. We 
consider three closest to the surface Cr monolayers with four 
types of Cr ions: site A (filled circle) and site B (hatched cir- 
cle) from the surface layer, site 2 (striped circle) and site 3 
(empty circle) from second and third closest to the surface Cr 
monolayers, respectively. Nearest-neighbor exchange param- 
eters between different types of Cr ions are denoted by thick 
gray lines. The red arrows show direction of local magnetic 
moment in the bulk-like AFM order. 

The parameter Ef M represents the surface energy in 
the paramagnetic state with no spin correlations. These 
energies are listed in Table HI] along with the configu- 
rational interaction parameters fitted to them. Again 
the three-body term is essential for obtaining a good fit. 
Without this term much larger CV score is obtained, and 
the take-one-out prediction for model B is about 43 meV 
off. As seen, the paramagnetic energies and the inter- 
action parameters differ little from their AFM values, 



TABLE III. Fitted parameters of Eq. © and misfits A (all in 
meV units) for surface models A, B, and A2B. The subscripts 
of exchange fields Hi and pair parameters Jij refer to the cor- 
responding Cr sites near the surface (see text). Last column: 
corresponding values in bulk Cr203 using values from Ref . I22I. 





A 


B 


A 2 B 


Bulk 


H A 


0.6 




0.9 


Jl = -2.2 


H b 




74.9 


69.1 


4 = "2.2 


H 2 


-37.5 


-16.8 


-29.1 


V b + 3Jl = -5.7 


H 3 


-1.4 


4.9 


-0.2 


3j£ + 3 Jl + Jl = 10.3 


JA2 


4.7 




6.4 


Jl = 2.1 


JA3 


11.1 




10.7 


4 = 3.0 


Jb2 




12.1 


4.8 


j| = 2.1 


Jbz 




5.4 


3.5 


Jl = 3.0 


J23 


-9.1 


0.6 


-8.6 


jl = -11.1 


A 


0.5 


0.1 


10~ 4 





indicating that magnctostructural coupling for the un- 
strained surface is weak. However, for the strained sur- 
face the parameters of the configurational Hamiltonian 
depend much stronger on the magnetic state, indicat- 
ing substantial magnetostructural coupling. Comparison 
of different columns of Table |H] shows that the effect of 
magnetic disorder for the strained surface is qualitatively 
opposite to that of the tensile strain. 

D. Ground-state phase diagram 

The search for the likely ground states of our model was 
performed by a direct enumeration of all configurations 
for unit cell sizes up to 6 x 6. The resulting ground-state 
phase diagram is shown in Figure [4] The Hamiltonian 
fitted to the surface energies of AFM unstrained super- 
cells (first column of Table H]) lies deep within the region 
where A2B (V3~ x \/3) is the ground state. Due to the 
long-range character of electrostatic interactions it is pos- 
sible that ground states with cell size larger than 6x6 
may appears in certain regions of the parameter space. 
However, such complicated orderings would be easily de- 
stroyed by thermal fluctuations. We therefore assume 
that such orderings, even if present, are irrelevant for the 
thermodynamic properties at temperatures where equili- 
bration is kinetically achievable; see discussion on kinetic 
energy barriers below. 



IV. CONFIGURATIONAL 
THERMODYNAMICS 

First we have studied the thermodynamics of our 
model within the mean- field approximation (MFA). We 
considered ordered structures A (1 x 1), AB (1x2), A2B 
(73x^3), and A 3 B 2 (y/3x ^7), which appear in the re- 
gion of the parameter space relevant for Cr203. We did 
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FIG. 4. Ground state phase diagram for the configurational 
Hamiltonian. The ground state configurations are shown in 
Figure [2] Red (green) circles and squares denote the values 
of parameters of the Hamiltonian fitted to the ground state 
(paramagnetic) surface energies for unstrained and strained 
Cr203 (0001) surface, respectively. The dashed line denote a 
strain path and the open black circles denote the parameters 
for intermediate strains. 



not include complicated orderings like A7B5 or A5B4, 
because, as noted above, they are expected to appear 
only at very low temperatures. The free energy of each 
phase was calculated in MFA; the equilibrium phase at a 
given temperature is the one with the lowest free energy. 
(Since the concentration of sites B is not conserved, the 
equilibrium phase is always single-phase.) 

The MFA results are shown in Fig. [5] (panels on the 
right-hand side). At small V/h or X/h, where the ground 
state is A, the surface never orders and remains A-type at 
all temperatures.-SA Where A2B is the ground state there 
is a continuous order-disorder transition from A2B-type 
to A-type. As the magnitude of V/h or X/h is increased, 
the critical temperature (in units of h) increases. This 
trend continues even when in the region where A3B2 
is the ground state. In this region, as temperature in- 
creases from zero, the A3B2-type structure undergoes a 
first-order transition to A2B-type, which then transforms 
to A-type at higher temperatures.— 

In the region where AB is the ground state, the sit- 
uation depends on X/h. For small X/h the AB-type 
structure undergoes a series of first-order transitions to 
AsB2-type and then to A2B-type; the latter then further 
transform to A-type. At larger X/h or V/h the A3B2- 
type ordering disappears and AB transforms directly to 
A 2 B. 

Due to strong geometric frustration, the MFA calcu- 
lations are unreliable and only provide a reference for 
comparison with more accurate calculations. We per- 
formed Monte Carlo (MC) simulations on L x L triangu- 
lar lattices with periodic boundary conditions. We usu- 
ally used L = 30 as this size is commensurate with all 
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FIG. 5. Temperature phase diagram for the configurational 
Hamiltonian obtained by MFA (right) and MC (left). Solid 
blue, dashed red, dash-dotted orange, and dotted green lines 
denote transition temperature to A2B (v3 x A3B2, 
A7B5, and AB orderings, respectively. In the MC case these 
lines are guides to the eye. Low-temperature MFA solutions 
are not shown in the patterned regions, because the corre- 
sponding ground states were not considered in the calcula- 
tions. 



the relevant orderings. In the loop over the lattice sites, 
a new state with the changed occupation number is tried 
and accepted or rejected using the Metropolis algorithm. 
The evaluation of the energy difference involves an ex- 
pensive calculation of the long-range interaction part (|5|) . 
For the two-body term this can be done using Fourier 
transforms. On the other hand, the direct calculation 
of the three-body term in Fourier space would be very 
expensive, because it requires a double summation over 
q. Instead, we calculate the Fourier transform of Ri as 
Rq = JqTiq, transform it back to real space using a fast 
Fourier transform (FFT) technique, and then calculate 
the three-body term in real space. We thus replace one 
sum over q by an FFT, which significantly reduces the 
computational cost for a large L. With this procedure 
the calculations could be performed for lattices with up 
to L = 36 using a few million (a few hundred thousand) 
MC steps per site for accumulating averages (for equili- 
bration). These restrictions were not always sufficient to 
obtain quantitatively accurate results (see below), but a 
qualitative understanding of the phase diagram could be 
achieved. 

The ordering type was identified by analyzing the 
structure factor I(q) = |n(q)| 2 , where n(q) is the Fourier 
transform of rij. All phase transitions between different 
ordered phases that we found are required by symmetry 
to be first-order. The order of order-disorder transitions 
was determined by analyzing the scaling behavior of the 
fourth-order energy cumulanti^ The transition tempera- 
tures were found from the peaks of the heat capacity for 
first-order transitions and from the finite-size scaling be- 
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havior of the fourth-order cumulant of the corresponding 
order parameter— for continuous transitions. For scaling 
analysis we used lattices with L = 30, 33, 36 for A 2 B-type 
ordering and L = 30, 32, 34, 36 for AB-type ordering. 

The results of MC simulations are shown in Fig. [5] 
(panels on the left-hand side). Similarly to MFA, for 
small values of V/h and X/h where A 2 B is the ground 
state an order-disorder transition to the A-type phase 
is observed. Our procedure identifies this transition as 
being everywhere continuous (second-order), except per- 
haps for small values of V/h, where the results suggest 
the proximity of a first-order transition. Note that the 
existence of a tricritical point was reported for the phase 
diagram of a related 2D hexagonal Ising model with AFM 
nearest neighbor and FM second-neighbor interactions. 35 

As expected due to strong geometric frustration, the 
critical temperature of the A2B-type ordering transition 
is strongly suppressed compared to MFA. For the param- 
eters corresponding to magnetically ordered unstrained 
Cr203 surface we found Tc = 165 ± 5i^T in MC compared 
to 600 K in MFA. If V/h or X/h are increased, initially 
Tc also increases due to the stabilization of the A2B- 
type structure relative to A-type. However, the increase 
of V/h or X/h also leads to stronger frustration, which 
tends to decrease Tc ■ This competition results in a max- 
imum of Tc as a function of these parameters. Note that 
the latter effect is absent in MFA (which is insensitive 
to frustration), which thereby completely fails for large 
V/h or X/h, incorrectly predicting that Tc should keep 
increasing. 

At low temperatures the system becomes difficult to 
equilibrate, and the equilibration time increases with in- 
creasing V/h or X/h. This is likely associated with in- 
creased frustration. In particular, we were unable to 
equilibrate the system at low temperatures in the pa- 
rameter range where the ground state is different from 
A and A2B, and the system remained in the initially 
chosen ordering state. In this case we chose the ini- 
tial state to be the ground state structure for the given 
set of parameters^ In particular, we performed MC 
simulations by increasing the temperature starting from 
the A 3 B 2 , A7B5, and AB ground states. Usually these 
structures underwent a first-order transition to A2B-type 
phase, which transforms to A-type under further heating 
(as discussed above). These two transitions are often 
very close to each other. Unfortunately, the temper- 
atures of both first-order and second-order transitions 
could usually be determined only with fairly large error 
bars. For first-order transitions these error bars are due 
to hysteretic behavior, while for the second-order tran- 
sition they result from strong fluctuations and limited 
averaging time. Often the error bars for these two tran- 
sitions overlapped, indicating that the appearance of the 
intermediate A 2 B phase might be spurious, and that the 
ground state structure may in reality transform directly 
to A-type. This is exactly what happens for large values 
of V/h and X/h when AB is the ground state. 



V. COMPARISON WITH EXPERIMENT 

Using the parameters fitted to the AFM surface ener- 
gies for Cr 2 03 surface the temperature dependence of the 
fraction of surface Cr ions occupying sites B was found 
from MC simulations, see Fig. [5] At low temperatures, 
when the surface is A2B-type, the concentration is close 
to the ideal value of 1/3 for the ground-state A2B struc- 
ture. As the temperature increases there is an order- 
disorder transition at Tc ~ 165K. In a temperature 
region around the transition the B-site fraction increases 
to about 40% and then stays approximately constant for 
temperatures well above room temperature. This result 
is in reasonable agreement with the room-temperature 
fraction of - 33% found from SXRD,1 

We found that the heat capacity has a broad shoulder 
above the critical temperature, indicating the persistence 
of strong short-range order well above room temperature. 
This is a direct consequence of geometric frustration^ 

The 1 X 1 — > -\/3 x y/3 ordering transition found above 
can be identified with the high-temperature phase tran- 
sition observed in LEED, 5 T c ~ 165 K produced by MC 
simulations agrees with the observed 5 - Tc ~ 150 K. How- 
ever, a second phase transition back to 1 x 1 at about 
100 K was also observed in these LEED measurements. 5 
This transition does not appear in our calculations for the 
unstrained surface. It was suggested 5 -^ that this second 
transition at 100 K may be induced by magnetostruc- 
tural coupling. As explained above, our calculations do 
not support this hypothesis for an unstrained surface. 
Our theory predicts that the surface of an unstrained 
Cr203 crystal should undergo only one 1 x 1 -> \/3 x vo 
ordering transition. 

For the analysis of the trend introduced by the strain 
we consider a continuous path in the parameter space 
assuming that all the parameters change linearly with 
strain, interpolating between the AFM surface energies 
found for 0% and 1.5% strain. As shown in Fig. @] the 
strain changes the ground state ordering from A2B to 
AB, passing through A 3 B 2 and A7B5 in between. The 
effect on structural thermodynamics is illustrated in Fig. 
[5J where the temperatures of different phase transitions 
are shown as a function of strain. The A2B-type ordering 
temperature is decreased by strain. At a certain value of 
strain the ground state changes to A 3 B 2 . Under heating 
this structure transforms to A2B, which then disorders 
at a higher temperature. For larger strains, however, the 
A2B phase disappears, and A3B2 transforms directly to 
A (1 x 1). As the strain further increases, the A7B5 phase 
is expected to appear at low temperatures (not shown in 
Fig. [5]) , and at yet a larger strain the AB phase sets in. 

The phase transitions at low temperatures may be un- 
observable for kinetic reasons. We have calculated the 
activation energy Eb for the jumping of a Cr ion from 
site A to B. Smooth barrier profiles were obtained with 
Eb equal to 0.4 eV and 0.3 eV for free and 1.5% strained 
surfaces, respectively. The frequency of thermally acti- 
vated jumps between sites A and B can be then estimated 
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FIG. 6. MC transition temperatures for different orderings 
as a function of tensile strain. Blue circles, red square and 
green rhombi denote transition temperatures below which 
A 2 B (a/3 x A 3 B 2 (a/3 x a/7), and AB (1 x 2), orderings 
respectively set in. Lines connecting the points are guides to 
the eye. 



the phase transitions are very sensitive to the growth 
conditions. This indirectly agrees with the fact that no 
phase transitions were observed for a thicker Cr2C>3 film. 9 
In this case the parameters of the Hamiltonian at low 
temperatures may correspond to a larger strain, which 
keeps the system disordered at all temperatures. 

The main drawback of the proposed mechanism of the 
reentrant phase transition is that magnetic disorder is 
assumed to influence the structural energetics at tem- 
peratures that are significantly below the Neel tempera- 
ture. Note, however, that the effect of magnetic disorder 
is to a notable extent mediated by the reduction of the 
parameter X, which describes the relaxation stiffness of 
site A relative to the electrostatic forces. But since the 
exchange coupling of site A to the bulk is quite weak (see 
Table UTT]) . the magnetic disorder affects this site already 
at low temperatures (see Section l"VIIj) . This factor gives 
some support to the proposed mechanism. For a more 
detailed consideration the magnetic degrees of freedom 
would have to be included in the Hamiltonian. We did 
not attempt this due to the limited amount of experimen- 
tal information on the atomic structure of the surface. 



as 7 ~ 7oe k s T where 70 is the attempt frequency on the 
order of a typical phonon frequency ~ 10 13 s _1 . (Or per- 
haps an order of magnitude smaller for thermal phonons 
at low T.) It follows that at room temperature the typi- 
cal hopping time is of the order of 10~ 8 s. On the other 
hand, the blocking temperature below which the kinetics 
is frozen is about 100 K. Therefore, the equilibrium phase 
transformations predicted for temperatures notably be- 
low 100 K are unobservable, and the system is expected 
to be trapped in the structural state corresponding to 
equilibrium near the blocking temperature. 

Our results support the hypothesis^ that magne- 
tostructural coupling plays an important role in the ori- 
gin of the two phase transitions observed in LEED for 
a thin strained film.™ The following picture can be sug- 
gested. At low temperatures the parameters of the con- 
figurational Hamiltonian correspond to the AFM-ordered 
strained surface. As seen in Fig. [6l in this case the 
equilibrium state near the blocking temperature is dis- 
ordered and has a 1 x 1 symmetry. Higher temperatures 
introduce partial spin disorder which, as discussed above, 
changes the parameters of the Hamiltonian similarly to 
a decrease of strain. This leads to the enhancement of 
the A2B-type ordering temperature (Fig. which be- 
comes higher than the blocking temperature and then 
overtakes the temperature of the system. In this picture 
this point corresponds to the low-temperature transition 
observed in LEED. As the temperature further increases, 
the system passes through the conventional disordering 
transition. 

The above scenario requires the surface to be under 
an exactly right amount of strain, and it implies that 



VI. SURFACE ELECTRONIC STRUCTURE 

In this section we discuss the electronic structure of 
the Cr203 (0001) surface. Fig. [7] shows partial densities 
of states (DOS) for A and B surface Cr ions for A, B, 
and A2B surface models. For comparison we include the 
partial DOS for the Cr ion in the middle of the A2B 
slab, which is similar to bulk CT2O3.— The DOS plots 
for different supercells are aligned using the semicore 2s 
states for bulk-like oxygen ions in the middle of the slabs. 

For model A2B there are two A-site Cr ions in the sur- 
face supercell, which are denoted as Ai and A2. The 
A2B ordering makes these sites inequivalent due to the 
directional character of the bonds, which was discussed 
in Section IlIIBI The partial DOS for Ai and A2 sites are 
similar, except for a shift of about 0.3 eV. This electro- 
static shift is due to the fact that site Ai is further away 
from the O sites in the subsurface layer than A2. 

The partial DOS for sites A in 1 x 1 and a/3 x a/3 sur- 
face supercells are qualitatively similar up to a moderate 
upward shift in the a/3 x a/3 supercell. The same can be 
said for site B, but the shift is in the opposite direction. 
The partial DOS for site A shows surface states in the 
bulk band gap both close to the valence band maximum 
and to the conduction band minimum. Site B introduces 
surface states originating from the valence band but ex- 
tending deeper into the bulk band gap. In both cases 
there is strong hybridization with the subsurface O ions. 
Partial DOS for deeper layers (not shown) shows that the 
surface states decay within 3-4 Cr monolayers from the 
surface. 
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FIG. 7. Spin resolved densities of states (DOSs) for A and B 
surface Cr ions and bulk-like Cr ion in the middle of the slab 
calculated for surface models: A (1 x 1), B (1 x 1), and A2B 
(v3 x y/3). Two nonequivalent A surface Cr ions for the A2B 
(v3 x y/3) surface model (see the text) are denoted by Ai 
and A2. Majority and minority DOSs are plotted on positive 
and negative y axis, respectively. Energy zero is set to the 
valence band maximum in the A2B (V3 x \/3) surface model. 
DOSs obtained for different slabs are aligned by semicore O 
2s states for bulk-like oxygen ions in the middle of the slabs. 
The green dashed vertical lines denote the bulk band gap. 
The DOS within the bulk band gap comes from the surface 
states. 



VII. SURFACE MAGNETISM 



cantly from the bulk couplings, which is a result of large 
ionic relaxations near the polar surface. (The parameters 
calculated as if the bulk exchange parameters^ do not 
change near the surface are listed in the last column of 
Table EH) 

To calculate the temperature dependence of magneti- 
zations for surface sites A and B, we used the mean-field 
approximation applied to the quantum spin-3/2 version 
of the Heisenberg model ©■ We considered the A 2 B sur- 
face model, since it is predicted to be the ground state for 
the unstrained Cr203 surface. Since the magnetostruc- 
tural coupling is weak, we expect that the surface site 
magnetizations are largely independent on the surface 
structure. We assumed that the exchange fields in © 
are proportional to the bulk mean-field sublattice mag- 
netization normalized to the experimental Neel tempera- 
ture. The resulting MFA surface site magnetizations are 
shown in Fig. [3] 
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In Section [V] we have seen that magnetic disorder may 
affect the surface phase transitions through a peculiar 
magnetostructural effect. On the other hand, the equilib- 
rium magnetization of the Cr203 (0001) surface enables 
interesting spintronic applications^ - — For these reasons 
it is interesting to consider the magnetic properties of the 
Cr203 (0001) surface at finite temperatures. 

The surface Heisenberg Hamiltonians were obtained in 
Section UTTl for different surface models (see Table Hill for 
the parameters for models A, B, and A2B). A common 
feature for all surface models is very strong exchange cou- 
pling of site B and weak coupling of site A to the bulk. 
This is expected, because all of the four bulk-like nearest 
and next-nearest neighbors of site A are absent; in spite 
of its large vertical relaxation, the remaining couplings 
do not compensate for this. The corresponding parame- 
ters for models A and A2B are quite similar. Although 
there are some differences for models B and A2B, con- 
figurations close to model B are statistically rare due to 
the fact that the equilibrium concentration of sites B is 
approximately 1/3. The fitted parameters differ signifi- 



FIG. 8. The temperature dependence of magnetizations of 
surface sites M(T) for the A2B (v3 x v3) surface model. 
Solid blue and red lines denote magnetization of site A and 
B, respectively. Dotted black line denote bulk sublattice mag- 
netization. Black circles show MC results for the temperature 
dependence of the concentration of surface Cr ions occupying 
site B. The solid black lines is the best fit to MC data. 

Since site B is strongly exchange-coupled to the bulk 
(Table IIIip . its magnetization largely follows the bulk 
sublattice magnetization. On the other hand, site A is 
weakly coupled to the bulk. As a result, its magnetiza- 
tion is substantially reduced and exhibits an inflection 
point. This inflection could be observed in the tempera- 
ture dependence of surface magnetic response in any ex- 
periment sensitive to the surface magnetism. Note that 
at 100 K site A is predicted to have already lost about 
60% of its magnetization at T = 0. As mentioned above 
in Section [V] this behavior lends some support to the 
magnetostructural coupling mechanism of the reentrant 
structural phase transition observed in Ref. [H. 



11 



VIII. CONCLUSIONS 

Based on first-principles total energy calculations and 
Monte Carlo simulations, we proposed a detailed micro- 
scopic model explaining the mechanisms of phase transi- 
tions at the stoichiometric Cr203 (0001) surface. Partial 
occupation of two surface Cr sites gives rise to compli- 
cated thermodynamic properties. Interaction is domi- 
nated by electrostatic forces, which promote ordering, 
and contains a smaller but still important contribution 
from atomic relaxations. The ground state is ordered 
with a v3 x y/3 unit cell; it undergoes a continuous 
order-disorder transition at Tq ~ 165 K. Tensile epi- 
taxial strain has a strong effect on the surface ener- 
getics, enhancing frustration, introducing new ground 
states and additional phase transitions. Magnetostruc- 
tural coupling also plays an important role in the struc- 
tural thermodynamics of the strained surface. Based on 
these results, we proposed an explanation of the reentrant 
1 x 1 -+ Vo x \/3 — > 1 x 1 phase transitions observed ex- 
perimentally on thin Cr2 03 (0001) films grown on Cr£ 



TABLE IV. Surface interlayer distances in % of the bulk in- 
terlayer distances for A (1 x 1) and B (1 x 1) surface models. 
The bulk interlayer Cr-O and Cr-Cr distances are 0.94 A and 
0.39 A, respectively^. Here A(n) denote nth atomic layer 
from the surface which has ions of type A. Our results are 
compared with existing literature. Here HF and MD denote 
Hartree-Fock and Molecular Dynamics methods, respectively. 
The experimental data (Exp) were obtained using LEED. 





A (1 x 1) 


B (1 x 1) 


LSDA+U 


GGA+IT 


HF b 


MD C 


Exp^ 


LSDA+U 


Cr(l)-0(2) 


-56.4 


-60 


-50 


-58 


-38 


-179.7 


0(2)-Cr(3) 


+7.3 


+ 12 


+3.3 





-21 


-9.2 


Cr(3)-Cr(4) 


-41.4 


-44 





-36 


-25 


-37.9 


Cr(4)-0(5) 


+10.8 


+9.2 





+17 


+ 11 


+18.5 


0(5)-Cr(6) 


+0.8 










+16.0 


Cr(6)-Cr(7) 


-2.4 










-45.4 


Cr(7)-0(8) 


+0.7 
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Appendix A: Surface relaxations 

Here we include the data on the atomic relaxations at 
the Cr203 (0001) surface and provide a justification for 
model ©-((U). Table ITVl lists the interlayer distances for 
A (1 x 1) and B (1 x 1) surface models. Strong inward 
relaxations are observed, as expected for a polar surface. 
For model A (1 x 1) the relaxations extend up to the fifth 
atomic layer, while for model B (1 x 1) they propagate 
much further, because the occupation of the interstitial 
site B introduces a stronger disturbance. Interlayer re- 
laxations for model A (1 x 1) are in reasonable agreement 
with other theoretical calculations & 12 ' 16 Although there 
are notable deviations from the LEED data^ we need to 
remember that the latter correspond to the actual sur- 
face termination but were fitted assuming the A (1 x 1) 
model. 

The upper panel of Fig.[9]shows the vertical coordinate 
of surface Cr ions occupying site A for different surface 
models as a function of Ri defined after Eq. (U). This co- 
ordinate is referenced with respect to that of the Cr ions 
occupying sites B averaged over different surface mod- 
els. (The subsurface O layer was used as an anchor for 
measuring the z coordinate.) For surface supercells with 
two inequivalent A-site Cr ions their vertical coordinates 



were similar, and we used their average. One can see that 
the calculated data agree very well with Eq. (Q} for both 
unstrained and strained surfaces. 



AAB(1x3) AAB(V3xV3) AB ABB (-/3xV3) ABB (1x3) 




FIG. 9. Upper panel: Vertical coordinate 2 of a Cr ion at 
site A for different surface models as a function of Ri defined 
after Eq. Q. Circles (squares) correspond to the unstrained 
(strained) surface. Solid lines are linear fits to the data. Lower 
panel: Basa function of R (see Eq. (|B2j) l for different surface 
models. Red (light) symbols correspond to the ground state, 
and blue (dark) symbols to the paramagnetic state. Circles 
(squares): data for unstrained (strained) surface. Solid lines 
are fits to a quadratic function with zero constant term, see 
Eq. (|B2|) . From bottom to top, the curves are shifted upward 
by 0, 0.01, 0.01 and 0.03 eV, respectively. 
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Appendix B: Quality of the fit 

Here we demonstrate the quality of the fit of ab initio 
energies to the configurational Hamiltonian and explain 
the importance of the three-body term in Eq. ([5]). 

Note that for all surface models for which ab initio en- 
ergies were calculated, the A sites are equivalent (ignor- 
ing the directionality of the bonds on the actual surface, 
see Section fill Bp . In this case, the surface energy from 
Eq. ([5]) can be rewritten as 

V o 
E = he- — (1 - c)(R + aR 2 ) + const (Bl) 

where R is the value of Ri for the A sites. Setting h — 
Eb — Ea, where Ea and Eb are the surface energies for 
models A (1 x 1) and B (1 x 1), we can define 

s ^E-E A -c(E B -E A)= _V V 

1 - c 2 2 v ' 



In the lower panel of Fig. [9] we plotted E as a function 
of R using the ab initio energies for all considered sur- 
face models. We included the data for both strained and 
unstrained surfaces using both ground state and para- 
magnetic surface energies. The resulting plots are very 
well fitted by the quadratic function with a zero constant 
term, demonstrating the high fidelity of the fit. 



The value of the parameter a extracted from the fit 
ranges from 0.04 to 0.05, which, as expected, is a small 
number. However, the relative importance of the three- 
body term compared to the two-body term is aR. Since 
for the considered surface models R varies between 4 and 
9, the relative importance of the three-body term is sub- 
stantial and reaches 50%. In the diagrammatic cluster- 
expansion language one can say that although a is small, 
the number of corresponding diagrams is large. 
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